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ABSTRACT 


One-  and  two-space  dimensional  finite-difference  schemes  for  the  Lagrangian 
numerical  solution  of  problems  in  the  motion  of  solids,  including  material 
strength,  are  presented.  IVo-dimensional  rectangular  cartesian  or  cylindri- 
cally  symmetric  problems  may  be  handled.  Results  of  sample  calculations  are 
appended  to  illustrate  the  effect  of  material  strength. 
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SECTION  I 


INTRODUCTION 


Solutions  of  problems  in  one-  and  two-dimensional 
(rectangular  cartesian  or  c> lindrically  symmetric)  fluid  flov; 
have  been  obtained  by  numerical  integration  using  high  speed 
digital  computers.  Such  methods  have  been  used  widely  to  solve 
problems  of  dynamic  deformation  of  solids  due  to  high  speed  im¬ 
pact  or  due  to  detonation  of  high  explosives.  The  assumption 
implicit  in  such  an  approach  is  that  material  strength  may  be 
neglected  in  comparison  with  the  high  pressures  attending  such 
motions. 

Experimental  evidence  on  craters  resulting  from  very 
high  velocity  particle  impact J  and  on  the  response  of  plates 
impacting  at  high  velocity^ indicates  that  material  strength 
does  influence  the  behaviour,  even  at  high  pressures  and  strain 
rates.  Aii  attempt  has  been  made  to  incorporate  material  strength 
in  a  two-dimensional  Lagrangian  numerical  integration  scheme. 
During  the  course  of  this  study,  two  similar  studies  were  re¬ 
ported.  The  present  work  is  essentially  similar  to  that  of 
Wilkins  and  Giroux^in  most  essential  features.  The  scheme  of 
Maenchen  and  Sack^is  based  on  a  different  method  of  obtaining 
finite  difference  analogs. 

The  relevant  differential  equations  have  been 
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collectedp  These  are  developed  into  finite  difference  equa¬ 
tions  by  the  application  of  Green's  Transformation.  Results 
of  a  few  sample  calculations  are  included  to  demonstrate  some 
of  the  qualitative  effects  of  material  strength. 
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SECTION  II 

DIFFERENTIAL  EQUATIONS 


2,1  Tensor  Equations 

The  governing  equations  are  the  differential  equations 

following  from  the  principles  of  conservation  of  mass,  momentum 
S 

and  energy.  From  the  principle  of  conservation  of  mass  follows 
the  continuity  equation 


P 


2.1 


where  p  is  the  density  of  the  deformed  volume  element  of  ma¬ 
terial  <sl  V  ,  and  the  initial  density  of  the  undeformed 

volume  element  of  material  d  \/  .  The  differential  form  of  the 
continuity  equation  is  also  useful. 


f 


2.2 


the  superimposed  dot  denotes  the  material  derivative,  and 
dj  is  the  stretching  tensor  defined  by  d  i,j)  • 

t*'"  =  3C '■  is  the  velocity  and  the  comma  denotes  the  covar¬ 

iant  derivative. 

The  equation  of  motion  follows  from  the  principle  of 
conservation  of  momentum 
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f 


2.3 


t 


£<••*3 


0 


where  o'**  U  **  is  the  acceleration,  ^  the  extrinsic 

(e.g.  gravity)  force  and  t  is  the  (symmetric)  stress  tensor. 

The  principle  of  energy  conservation  and  the  definition 
of  internal  energy  lead  to  the  energy  equation 

f  ^  *  f  Q.  2.4 

where  c  is  the  internal  energy  per  unit  mass,  h  is  the 
heat  flux  vector  and  d  the  (chemical  or  other)  energy  supply 
per  unit  mass.  It  is  more  convenient  to  resolve  the  first  term 
on  the  right,  representing  the  stress  work,  into  spherical  and 
deviatoric  parts.  Writing 

ti  -  i  ti  s:  *  "t:  2.5 

where  the  superscripted  d  denotes  the  deviator,  and  similarly 
for  d  ,  we  have 
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2.6 


t 


k  rn 


knt 


Writing  P  =  ~  ^3  where  p  is  the  pressure,  and  using  the 

differential  form  of  the  continuity  equation,  the  energy  equa¬ 
tion  takes  the  form 


f 


i 


■f 


k<n 


a 


2.7 


These  equations  are  supplemented  by  the  constitutive 
equation,  describing  the  material  behaviour,  and  the  boundary 
conditions,  which  are  described  later. 


2.2  Two-Dimensional  Equations  in  Physical  Components 

The  tensor  equations  are  expanded  in  terms  of 
physical  components  for  rectangular  and  cylindrical  polar  co¬ 
ordinates  for  two  independent  space  variables. 

In  two  independent  space  variables,  we  consider  rec¬ 
tangular  symmetry,  with  no  variation  in  quantities  along  the  y 
axis,  and  cylindrical  symmetry  with  no  variation  in  quantities 
in  the  tangential,  or  fi  direction.  The  physical  equations 
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for  both  cases  may  be  written  do\m  simultaneously  by  v^n^iting  x- 
and  tj  in  place  of  r  and  9  for  the  cylindrical  polar  case, 
and  defining  o<  =  /  for  rectangular,  cx  =  2  for  cylindri¬ 
cal  symmetry. 


The  equations  of  motion  become 


? 


= 


at 


%  X 


-  t 


3  >c 


f 


a 


=  pf 


rt_ 

3  X 


xz 


Z2 


t 

(oc-l)  -i_ 

^  ''  5C 


yy 


)  2.8 


The  velocity  gradient  Uij  has  physical  components 


u 


«,b 


ax 


dx. 


/  ^  ^ 


0 


az. 

0 

az 


2.9 


where  a  ,  b  take  on  successive  values  x  ,  ^  ,  z  and  it  is  im¬ 
material  whether  indices  are  written  as  subscripts  or  superscripts 
since  the  coordinates  are  orthogonal. 


The  differential  form  of  the  continuity  equation  be¬ 
comes,  directly 


6 


3u 
3  X 


a—  j 


2.10 


In  order  to  expand  the  energy  equation  we  require  the 
physical  components  of  the  stretching  deviator 


d:  -  T  <^1  s 


t  k  I  C  k 

^  ^  J  JD  ^  m 


>  2.11 


Thus  we  have  physical  components 


(ab) 


hi"  i  £ 

3x  ^  -3  f 


0 


I  1  )  p 


and  the  energy 

equation 

becomes 

if 

• 

p 

dp 

I  a” 

where 

the  deviator  stres 

3  work 

i  i  <?  z  a  >c  / 


0 


3u  ^  ^  jL  £ 
32,  ^  3  p 


d: 


2.12 


2.13 


7 


-  z  ‘'r*  ‘U"  +  I 


^  •'r"  ^  •u*’' 


2.14 


Since  material  particles  will  be  followed  in  the 
computation,  the  velocities  and  accelerations  are  simply 


given  by 


u 


u 


H 


2.15 


»  z 

a  -  a 

2.3  One-Dimensional  Equation^  in  Physical  Componc  ;ts 


The  tensor  equations  are  expanded  in  terms  of  phys¬ 
ical  components  for  rectangular,  cylindrical,  and  spherical 
polar  coordinates  with  variation  only  in  the  X  or  radial  di¬ 
rection.  Proceeding  similarly  as  in  the  two-dimensional  case, 
further  defining  **<>=■  7  for  spherical  symmetry.  Then  the  equa¬ 
tion  of  motion  reduces  to 

yt  /)K  ^  t  * 

rf  ^  ^  ~  2. 


There  is  no  rotation,  and  the  above  stress  components 
are  principal  components,  which  is  expressed  by  use  of  a  single 
superscript.  Off-diagonal  stretching  components  also  disappear, 
and 
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■  0  ■,  0] 


0<  e  I 


r 

I  aK  > 


u  ‘  1 

Z  >  o] 


jn  r  .  ii’  u.*] 

^  *  I  av  >  5c  '  X  J 

The  energy  equation  reduces  to 

aK* 


o<  =  i 


o<  *  3 


do  /  ,  A.*  1 

J  ^  ;  y-o. 


where  the  deviator  stress  work  is 


=  iV  ‘U* 


J|Z  ^  ^  M.m 


2**^  12 


2.17 


2.18 


2.19 
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SECTION  III 


CONSTITUTIVE  EQUATION 

3ol  Compressible  Fluid 

We  first  discuss  the  constitutive  equation  of  a  fluid, 
which  cannot  support  a  shear  stress  while  at  rest.  The  total 
stress  tj  entering  the  conservation  equations  may  be  expressed 
in  terms  of  a  thermodynamic  (elastic)  pressure  gp  ,  a  dissi¬ 
pative  (viscous)  pressure  ,  and  a  viscous  stress  deviator 

‘j  "  -(^F  *  f)  ^  Vi  3.1 


The  viscous  stress  tensor  must  be  isotropic.  Taking  the  vis¬ 
cous  stress  proportional  to  the  stretchings,  we  have  the  follow 
ing  form 

^  3.2 

where  ^  and  yU  are  here  used  as  viscosity  coefficients. 
Then 


3.3 
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and  for  the  deviator 


<4 


3.4 


Using  Equations  3.2  and  3.3  this  reduces  to 


A 


3.5 


A  ,  i 

where  a  •  is  the  previously  defined  stretching  deviator. 

The  viscosity  coefficients  may  be  "real"  coefficients, 
representing  the  viscous  behaviour  of  the  material,  if  these  can 
be  evaluated.  More  usually  artificial  viscosity  coefficients  are 
used  with  specific  properties  to  aid  in  maintaining  stability  of 
"he  finite  difference  calculation  without  unduly  affecting  the 
solution.^  These  will  be  discussed  later. 


The  thermodynamic  pressure  is  related  to  the 
thermodynamic  state  (  £  ).  This  relation  will  be  taken  in 
the  form  of  power  series  in  the  compression  7  * 


eP  ^  f,  IP>  *  <?•  /x  (P) 


-  K  ^  {  I  ^  h.  ^  ^  .  ...  j  3.6 

/x  '  *  ^'2  *■  7*  - } 

A  similar  enuation  of  state  is  used  for  explosion  product  gases, 
but  with  different  forms  of  the  functions  f^  and  f.^. 
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3.2  Elastic  Perfectly  Plastic  Material 
3.2.1  Tensor  Equations 


For  a  compressible  material  which  is  able  to  support 
a  shear  stress  while  at  rest,  it  is  necessary  to  add  a  non- 
dissipative  (elastic)  stress  deviator  to  the  previous  expression 
for  the  total  stress. 


t 


J 


3.7 


The  thermodynamic  pressure  and  viscous  stresses  will  be  formu¬ 
lated  as  in  the  case  of  a  fluid.  Here  we  will  consider  the 
elastic  stress  deviator. 


Making  the  usual  assumptions  of  plasticity  theory, 
modified  to  include  large  compressions’  viz:- 

a)  The  stretching  can  be  expressed  as  the  sum  of  the 
elastic  and  plastic  stretchings. 


3.8 


b) 


elastic  and 


The  spherical  stretching  (dilatation)  is  entirely 
recoverable,  which  is  equivalent  to  assuming 


3.9 


c)  The  elastic  stress  deviator  is  limited  by  the  von 
Mises  yield  condition 


12 


3.10 


(p^S-jy^io 

where  3  ■  is  the  second  moment  of  the  stress 

deviator,  and  Y  is  a  material  constant.  It  is  easily  seen 
that  y  is  the  yield  stress  in  a  uniaxial  stress  tensile  test. 

d)  The  plastic  stretching  is  orthogonal  to  the  yield 
surface,  leading  to  the  flow  rule 


f  J 


3.11 


where  JT  is  an  undetermined  proportionality  constant.  For 
von  Mises  yield  criterion,  this  reduces  directly  to 


p  J 


3.12 


In  view  of  the  fact  that  we  have  assumed  ^  0  y  vje  may  also 

write  the  stretching  deviator  on  the  left  hand  side  of  this 
equation . 

e)  The  usual  assumption  of  a  linear  relation  between 
elastic  stress  rate  and  stretching  must  be  modified  to  allow 
for  finite  compression.  The  usual  assumption  of  an  isotropic 
linear  elastic  medium  is  expressed  in  Hooke's  Law 

tj  =  Ad;  SJ  ♦  2^  A]  3.1 


13 


v;hcrc  X  and  ^  are  now  used  as  elasLicily  (Lame)  consLanLs, 
and  the  superimposed  Lriangle  denotes  ooIccllvc  stress  rate. 
Resolvinj^  this  expression  into  spherical  and  dcviatoric  parts 
as  was  done  lor  the  viscous  stress  tensor, 


-  P  ‘  J 

't;  -  V 


3.14 


We  nov;  make  tlie  special  assumption  that  the  clastic  devi^itoric 
strain,  which  is  limited  in  effect  oy  the  yield  condition,  alv;ays 
remains  small.  We  are  then  concerned  only  v;ith  small  deviations 
from  a  hydrostatic  state,  and  :.iake  the  assumption  tliat  the  ma¬ 
terial  remains  isotropic.  We  thus  generalize  the  above  equa¬ 
tions  by  writing 


f  . 


2  y 


3. 15 


where 
Lalcen 
It  is 
ol  liq 


K  ,  the  bulk  modulus,  and  d  ,  the  shear  modulus,  arc 
t(2  iic  functions  of  the  thermodynamic  state  f  P,  ^ ' 
of  course  more  convenient  to  use  the  integrated  version 
nation  3.13a,  that  is  Equation  3.6,  for  the  spherical  part 
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The  equaLions  accessary  Lo  dcLerniine  Lhe  elasLic  stress 
deviaLoi'  ''ro;,!  clu'  sLrc'Lchin;’  deviaLor  are  LhereL'ore 


J 

f 

P 


J 


(deco!;ipt)s  i  L  ioa) 


2  G 

e  J 


(l' las  Lie  re  1  alien) 


3.  10 


E 


$ 


(;■  iL  ]  Cl  erili’rion) 


'’a' 

p  J 


(ilow  rule) 


In  the  next  two  sections  the  second  and  third 
equations  ot  the  above  set  will  be  expanded  into  physical 
components.  The  first  and  fourth  equations  pose  no  difficulty. 

3.2.2.  Tv;o-Dlriiensioaa  1  Equations  in  Physical  Coniponenls 

lixpaniling  Liu;  siconil  iuonc-nL  of  llic  siren'.:',  deviator  in 
the  Lv;o-dimcnsional  case,  the  nonzero  terms  are 


u 


3.17 


IS 


where  use  has  been  made  of  the  fact  that  the  trace  of  a 
deviator  vanishes. 

The  definition  of  objective  stress  rate  is 


r 

t 


i.r 


3.18 


where  uT;;  is  the  spin  tensor  defined  as  uT/.- =  Ur:  .-t  .  Refer- 
ring  to  the  physical  components  of  the  velocity  gradient  given 
previously,  the  only  nonzero  component  is 

au*  au* 


ur 


X  z 


®  i  I  az  a>c  / 


3.19 


so  that  the  physical  components  of  the  objective  stress  rate 
for  both  rectangular  and  cylindrical  coordinates  are 


« 


at“ 


z  u/"*  it 


X  z 


2  G 
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It  is  unnecessary  to  refer  to  principal  directions  in 
order  to  determine  the  stress  deviator.  If  principal  components 
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and  directions  are  needed,  for  example,  to  apply  a  fracture 
criterion,  these  can  be  found  very  simply.  We  note  that  the 
stress  deviator  matrix  is 


so  that  one  principal  direction  is  the  y  direction  with 
principal  component  and  the  other  two  principal 

components  are  found  by  diagonalizing  its  cofactor 


=  1  *  (yy  )  3.21 


where  a  -  1,2.  respectively. 

Components  of  the  unit  vectors  corresponding  to  these  principal 

values  are  then  given  by 
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These  components  are  direction  cosines  of  the  principal  direc¬ 
tions.  It  is  easily  seen  that  the  angles  between  the  principal 
directions  and  the  x  axis  are 
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where  the  two  values  of  ^  from  equation  3.21  are  used. 

3.2.3  One-Dimensional  Equations  in  Physical  Components 

In  one  dimension,  the  equations  are  somev;hat  simpler 
due  to  the  absence  of  rotation.  The  yield  condition  reduces  to 
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while  the  clastic  relation  becomes. 


H 


2  /d  * 


For  the  rectangular  (  (x  *  / 
cases,  ehe  symmetry  of  the  motion  leads 
cation.  Since 

d.Y  _  d.z. 


)  and  spherical 
to  considerable 
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and  Lhe  conciilion 
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the  yield  condition  iminediatcl)'  reduces  to 
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It  is  Lhcreforc  unnecessary  to  refer  to  t’ne  rest  of  equations 
3.16  when  the  material  is  at  yield.  When  the  material  is  elas- 
tic,  equation  3.25  suffices  to  determine  . 


For  the  cylindrical  case  (  ot  -  2.  )  such  simplification 
is  not  possible  It  is  thus  necessary  to  solve  the  v;hole  sys¬ 
tem  of  equations  3.16.  The  clastic  relation  in  physical  com¬ 
ponents  becomes 
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l!;  is  nnnocessary  to  use  the  second  relatluii  above.  using  the 

y 

property  that  the  trace  of  a  deviator  vanishes  to  eliminate 
from  the  yield  condition, 
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3.2,4  Geometrical  Representation  of  Stress 
If  the  stress  deviator  is  referred  to  principal  axes, 
which  in  the  one-dimensional  case  coincide  with  the  coordinate 
axes,  the  shear  components  vanish.  It  is  then  possible  to  plot 
the  stress  state  in  a  three  dimensional  rectangular  coordinate 
systemf  The  condition  that  the  trace  of  the  stress  deviator 
vanishes,  equation  3.27,  defines  a  plane,  called  the  7T  plane, 
with  a  normal  which  has  direction  cosines  ( y/r,  'y/r,  '4r)- 

The  yield  condition,  equation 
3.24  defines  a  sphere  of  radius 
JJy  ,  Thus  all  attainable 
stress  deviator  states  are  lim¬ 
ited  to  the  7T  plane  in  a  domain 
within  a  yield  circle  of  radius 


/py 


It  is  convenient 


^  to  rotate  the  coordinate  system 
so  that  the  yield  circle  is  in 
the  plane  of  the  paper,  and  to 
use  polar  coordinates  to  deter¬ 
mine  any  given  stress  state  P 
in  this  plane. 
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It  is  iinmediately  evident  from 
equation  3.24  that  the  radius 
vector  is  given  by  /F.  The  an¬ 
gle  between  the  radius  vector, 
and  the  line  formed  by  the  inter¬ 
section  of  the  TT  plane  and  the 
coordinate  plane  is 

given  by 

(-  M')  3.31 

where  2  ^  *  -  ‘‘t*' 

»  — - i - 

•'f  - 

is  Lode's  variableo  Using  equation  3.27,  this  reduces  simply  to 

-f  .  »rct<in  3.32 

'  ■‘t"  - 

It  is  clear  that  J[  ,  and  p  are  sufficient  to  define 

any  elastic  stress  state. 
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SECTION  IV 


STABILITY 


4.1  Artificial  Viscosity 

Discontinuities  or  shocks  may  develop  in  solutions  to 
the  motion  of  perfect  compressible  fluids  and  solids.  Such  dis 
continuities  lead  to  instabilities  in  the  finite  difference  cal 
culationv.  This  problem  is  avoided  by  rendering  the  solution 
smooth  and  continuous  everywhere  by  the  introduction  of  artifi¬ 
cial  viscosity,  so  formulated  that  solutions  are  only  affected 
in  areas  of  very  high  gradients,  i.e.,in  shock  zones.  Follow- 
ing  von  Neumann  and  Richtmyer  we  choose  a  bulk  viscosity  coef¬ 
ficient  dependent  on  the  dilatation  so  that  the  viscous  stress 
is  negligible  in  areas  of  moderate  gradient,  i.e., 

*  f  k 

in  equation  3.3,  where  is  a  constant  with  dimensions  of 

length.  In  some  problems  it  is  found  desirable  tc  include  a 
linear  viscous  coefficient.  Following  Landshoff  this  is  for¬ 
mulated  as 

(a  ♦  c  4.2 
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where  A  also  has  dimensions  of  length,  and  c  is  the  sonic  ve¬ 
locity  of  the  medium.  The  viscous  pressure  can  therefore  be 
written,  with  the  use  of  the  continuity  equation,  as 


.  1 


4.3 


Since  only  compression  shocks  are  possible,  (y  is  set  equal  to 
zero  when  p  k  O 

Experience  has  shown  that  equivoluminal  oscillations 
may  occur  in  two-dimensional  finite  difference  schemes,  and  a 
variety  of  artifices  have  been  adopted  to  strengthen  the  viscous 
stress  in  the  direction  of  maximum  velocity  gradient.  This  can 
be  done  elegantly  by  introducing  a  viscous  stress  deviator  in  a 
way  analogous  to  the  introduction  of  the  viscous  pressure  above. 


4.4 


Again,  the  viscous  stress  may  be  unnecessary  on  expansion, 
and  may  be  set  equal  to  zero  when  >  0  . 

There  is  no  difficulty  in  expanding  equation  4.4 
into  three  equations  in  physical  components  for  the  two-dimensiona 
case.  In  the  one-dimensional  case  it  is  sufficient  to  retain 


only  the  bulk  viscosity. 


4.2  Stability  Criterion 


The  conditional  stability  of  second  order  finite  dif¬ 
ference  equations  is  well  knownf  In  the  finite  difference 
scheme,  quantities  are  sampled  at  discrete  intervals  in  space 
(  Ax  ,  A  a:  )  and  time  t  ),  The  time  increment  cannot 

be  chosen  arbitrarily.  If  the  time  increment  is  too  large, 
disturbances  with  wave  lengths  of  the  order  of  the  mesh  size 
tend  to  grow  without  bound.  A  variational  analysis  of  the  one 
dimensional  finite  difference  equations  (Appendix  A)  in  which 
solutions  are  sought  in  the  form 
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etc.  suggests  that  a  sufficient  condition  for  stability  is 

"  (\¥Xf3^)c  r 

Deviator  viscous  stresses  have  been  omitted  in  this  development, 
since  they  are  in  general  smaller  than  the  spherical  viscous 
terms.  No  difficulty  has  been  encountered  on  this  account,  pre¬ 
sumably  because  equation  4.6  is  somewhat  more  stringent  than 
the  necessary  condition. 


A  complete  stability  analysis  has  not  been  carried  out 
for  the  two-dimensional  finite  difference  equations  (Appendix  B) . 
However,  an  analogous  form  of  stability  criterion  has  been  found 
useful.  Note  that  in  one  dimension  the  continuity  equation 
can  be  written  (equation  2.2) 
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Taking  as  a  measure  of  mesh  size,  where  f\  is  Che  area  of 


a  mesh,  (^quaCion  4.6  becomes 


*  j.  ^  ^ 

^  ~ — : — j — i  I  4.8 

Despite  the  fact  that  /T  is  not  a  good  approximation  to  the 
mesh  size  for  a  highly  distorted  mesh,  the  use  of  this  equation 
has  not,  so  far,  led  Co  any  stability  problems.  This  presum¬ 
ably  arises  from  the  fact  that  the  error  in  the  numerator 
partially  offsets  the  error  in  the  denominator. 
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SECTION  V 


FINITE-DIFFERENCE  ANALOGS 

The  differential  equations  given  in  the  previous 
section  are  to  be  integrated  numerically.  The  derivative 
terms  are  replaced  by  finite  difference  analogs  to  produce 
a  set  of  algebraic  equations  which  are  solved  in  a  stepwise 
manner  to  produce  the  desired  solution.  The  principal  dif¬ 
ficulty  concerns  the  choice  of  finite  difference  analogs 
to  the  partial  derivatives. 

In  one  space  dimension  centered  difference  ex¬ 
pressions  correct  to  second  order  will  be  used.  In  two 
space  dimensions,  a  variety  of  difference  analogs,  presum¬ 
ably  of  the  same  order  of  approximation,  have  appeared  in 
the  literature.  Some  of  these  have  been  compared  in  another 
report.***  Here  short  derivations  of  the  analogs  which  were 
found  to  be  most  suitable  are  given.  Application  of  these 
to  develop  the  complete  finite  difference  equations  in  one 
and  two  space  dimensions  is  postponed  to  the  appendices. 

In  one  dimension,  the  usual  approach  is  to  use 
Taylor's  expansion.  Consider  an  arbitrary  quantity 
varying  with  an  independent  variable  x  (which  may  represent 
time,  or  distance  in  one  space  dimension).  We  wish  to  find 
the  gradient  of  at  a  point  0  given  values  of  at  neigh¬ 
boring  points  1  and  2  finite  distances  from  0.  Applying 
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Taylor's  expansion, 


d  \y  I  /  \  “ 


X 


5.1 


,x 


5.2 


Solving  these  for  , 

if.  .J'if'  r  -  iff  f  (f  i-IfiT  ,  ji  j  ]-■■■ 

where  f  J  ~  ^  measure  of  the  assym- 

metry  of  the  mesh.  If  the  mesh  is  nearly  symmetric  the  second 
order  term  is  negligible,  and  if  the  mesh  size  is  small,  the 
higher  order  terms  are  negligible,  and 

if  .  Jtv.-  f,  5.3 

Xx  - 

It  is  clear  that  the  truncation  error  grows  as  the 
mesh  becomes  more  asymmetric. 

Precisely  the  same  reasoning  has  been  used  to  con¬ 
struct  a  finite  difference  analog  uo  partial  derivatives  in 

two  space  dimensions.  Consider  an  arbitrary  quantity  ^  vary- 

f 

ing  with  two  independent  vari¬ 
ables  X  and  z.  We  wish  to  find 
the  gradients  of  in  the  x  and 
directions  at  point  0  given 
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values  of  at  four  neighboring  points  1234  finite  dis¬ 
tances  from  0,  Applying  Taylor's  theorem  we  obtain  four  equa¬ 
tions  of  the  type 


^  ID 

y',  -  fo  ♦  ’■  ■*»>'  jr 


Kolsky^'suggested  solving  this  overdetermined  system  of  equations 
by  first  solving  for  (  </',  -yj.)  and  an.  -  K'*  )  giving  two 

equations 


where  the  remainder  terms  again  depend  on  the  assymmetry  of  the 
mesh  and  the  mesh  size.  Again,  providing  that  the  ass5nnmetry 
and  mesh  size  are  both  small,  the  remainder  terms  may  be  neg¬ 
lected  in  comparision  with  the  first  order  terms.  Solving  for 


t-he  gradients  in  this  case 

=  IJ  {(n.- 
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where 


It  is  seen  that  A  represents  the  area  of  the  quadrilateral 
1  2  3  4. 

It  is  interesting  to  note  that  an  identical  result 
may  be  obtained  in  another  way.  Green's  Transformation  in 
two  dimensions  may  be  written 

f  >i!,i  da  ^  f  a  ds  5.7 

5 

where  n-  is  the  unit  exterior  normal  to  the  surface  S  er^- 
closing  an  area  A.  In  component  form  this  becomes 


S  /? 


where  the  gradients  have  been  averaged  over  the  area  A.  Apply 
ing  these  relations  to  the  quadrilateral  1234  surrounding 
point  0,  we  find  that  the  average  gradients  may  be  expressed  a 

^  =  J  f  !/'rx  -^*J  } 

5.9 
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where  <i’r.  is  the  average  of  over  the  side  1  2  etc.  if  we 

write 

‘  i  [  IC,  ^ 

the  above  equations  5.9  immediately  reduce  to  the  previous 
equation  5.6. 

Both  forms  (equations  5.6  and  5.9)  will  be  found  use¬ 
ful  in  developing  the  finite  difference  equations  in  two  space 
dimensions.  It  is  clear  from  their  development,  however,  that 
the  truncation  error  due  to  neglect  of  the  higher  order  terms 
must  depend  on  the  asymm'^try  of  the  mesh.  The  truncation  error 
has  been  investigated  in  another  report,  and  has  't)oen  found  to 
become  comparable  to  the  terms  which  are  retained  even  for  mod¬ 
erate  distortions  of  the  mesh. 
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SECTION  VI 


RESULTS 

A  few  illustrative  calculations  are  |:)resented  in  this 
section.  While  the  results  are  suggestive  of  the  role  of  mate¬ 
rial  strength  in  plastic  v/ave  propagation,  considerable  work 
remains  to  be  done,  particularly  on  tlu-  reprc'sen La t i on  of  ::iaU-- 
rial  properties,  before  cjuant  i  tat ive  information  can  no  extracted. 

Results  using  the  one-d  Lmens  Lena  1  fiixiLe  differL>ncL‘ 
method  for  Llie  faco-on  impact  of  tWvi  aiuMiniim  plates  are  shown 
in  Figs.l  through  4.  The  configurations  were  chosen  to  cor- 

X 

respond  to  experimental  data  reported  by  Curren.  The  materi¬ 
al  constants  used  were 

/•  =  2.78  gm  cm ^  V  =  2.7b  kb 

=  764  kb  =  3.473  fc  =  286  kb 

=  1  014  ^  =  1.0  =  5.86 

kx  =  -.  2  36  =  1.0  =  8.44 

kj  =  -.513  =  1.0  fi  =  0.0 

Figure  1  shov.’s  print  plots  of  stress  profiles  for  a 
driver  plate  velocity  of  1.9  Ion  sec,  while  1-ig.  2  shows  corre¬ 
sponding  results  when  the  yield  stress  is  set  to  zero.  Figures 
3  and  4  show  initial  rear  surface  velocities  of  the  target  plate 
upon  reflection  of  the  stress  wave  as  a  function  of  target  plate 

thiclaiess  (in  terms  of  tlriver  plate  thickness),  for  driver  plate 
velocities  of  1.9  1cm  sec  and  1.2  Icn  sec  rc'spec Lively .  Also 

shown  are  experimental  data  of  Curran,  and  results  of  the 
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hydrodynamic  r.heory  of  Fowlcs.*^ 

The  computer  solutions  for  zero  strength  {  y  =  0  ) 

fall  well  below  Fowles' prediction.  This  may  be  ascribed  to  dif¬ 
ferences  in  the  fit  to  the  hugoniot,  and  the  fact  that  Fowles 
did  not  account  for  entropy  changes,  in  so  far  as  he  used  the 
hugoniot  for  the  expansion  instead  of  an  isentrope.  The  com¬ 
puter  solutions  including  material  strength  fall  well  belcw  the 
zero  strength  solutions,  i>ut  not  sufficiently  to  agree  v;ith  the 
experimental  data.  The  remaining  disitgreement  may  well  be  due 
tj  an  increase  in  '’ield  strength  with  compression  as  suggested 
bv  Curran. 

ir/lien  material  strength  is  included,  it  is  seen  from 
Fig.  1  that  an  clastic  release  wave  with  an  amplitude  of  about 
15  kb  moves  into  the  compressed  material  behind  the  shock  with 
the  local  elastic  wave  velocity.  The  existence  of  a  15  kb 
elastic  release  wave,  despite  the  fact  that  the  yield  stress  is 
taken  constant  at  2  76  kb,  can  be  explained  by  referring  to 
Fig.  5  which  shows  a  schematic  of  the  clastic  plastic  constitu¬ 
tive  relation  neglecting  hysteresis  due  to  the  entropy  change  in 
the  shock.  On  loading,  a  stress-strain  path  lying  a  distance 
yj  V  above  the  hydrostatic  curve  is  followed  (0  A  B  C  in 
Figo  5).  On  release  of  stress  from  state  C,  path  C  D  E  is  fol¬ 
lowed,  so  that  an  elastic  release  wave  of  amplitude  C  D  is  gen¬ 
erated.  The  amplitude  of  the  elastic  release  wave  is  determined 
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by  the  slope  of  segment  C  D,  relative  to  the  slope  of  the  hy¬ 
drostats.  The  slope  of  C  D  is  (  K  ^  )  while  tlie  slope 

of  the  hydrostat  is  (^).  Thus  the  relative  slope  depends  on 
the  shear  modulus  6  .  At  the  same  time  the  velocity  of  th.e 
elastic  release  wave  is  given  by 

Thus,  if  the  value  of  Ca  is  increased,  the  elastic  wave  amplitude 
is  decreased  but  its  velocity  is  increased.  Higher  attenuation 
results  from  a  larger  clastic  release  wave  amplitude,  but  also 
from  a  release  wave  of  higher  velocity.  Thus  the  two  effects 
offset  each  other  to  some  extent,  and  some  variation  in  (»  does 
not  materially  affect  the  results,  as  found  by  Jones/ 

Once  the  elastic  release  wave  has  reached  the  loading 
shock,  the  shock  amplitude  is  reduced.  As  the  reduced  amplitude 
shock  propagates,  the  material  is  now  loaded  to  a  lower  stress 
(state  B,  Fig.  5).  Wlicn  the  release  wave  reaches  this  material, 
a  new  elastic  release  wave  will  appear  (corresponding  to  segment 
B  E,  Fig.  5)  and  the  attenuation  process  will  be  repeated. 

Results,  using  the  one-dimensional  finite  difference 
method,  for  the  response  of  an  aluminum  sphere  18  cm  in  diameter 
containing  a  concentric  spherical  cavity  3.4  cm  in  diameter  fil¬ 
led  with  Pentolite  are  shown  in  Fig.  6  (not  all  mesh  points  arc 
plotted).  The  configuration  was  chosen  to  correspond  to  experi- 

ir 

mental  data  obtained  at  the  Ballistics  Research  Laboratory.  The 
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material  constants  for  the  aluminum  were  identical  to  those  used 
above,  while  the  material  constants  for  the  Pentolite  were  taken 
to  be 

=  1.714  gm/cm^  =  2.333  gm/cm^ 

r  =  1.77  p,  =  290.289  kb. 

=  2.77 

D  =  7.991  lon/'sec 

It  can  be  seen  that  the  results  for  zero  strength 
{  O  )  and  a  strength  Y  =  2.76  kb  do  not  differ  very  materially 
except  at  late  times.  The  difference  is  sufficient  however  to 
cause  a  drastic  difference  in  spall  behaviour.  The  elastic- 
plastic  case  (y=  2.76  kb)  showed  3  spalls  at  radii  of  7.55, 

5.70  and  3.90  cm.  respectively,  while  the  hydrodynamic  case 
(y»0)  shov;od  several  spalls,  the  outermost  occuring  at  a  radius 
of  3.55  cm.  Great  care  is  necessary  in  interpreting  spall  results 
from  finite  difference  calculations,  since  these  depend  to  some 
extent  on  the  choice  of  mesh  size  and  artificial  viscosity  co¬ 
efficients.  The  auovc  results  should  be  regarded  as  preliminary. 

In  Fig.  7  is  plotted  the  cavity  radius  vs.  time  for 
y  =  0  and  y=  2,73  kb.  The  cavity  radius  grows  somewhat  more 
slowly  when  material  strength  is  included.  The  cavity  radius 
at  20  microseconds  was  computed  to  be  2.68  cm,  which  might  be 
compared  with  the  final  cavity  radius  of  about  2,95  cm  found  in 
the  experiment.  It  is  clear  that  the  cavity  grov;s  for  times  in 
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excess  of  100  microseconds ^  which  is  lonj^  compared  to  the  time 
required  for  the  initial  shock  to  reach  the  outside  boundary  of 
the  aluminum  (about  15  microseconds) . 

Preliminary  results  using  the  two-dimensional  finite 
difference  method  for  the  end-on  impact  of  a  finite  length 
aluminum  cylinder  on  a  smooth  wall  (or  symmetry  plane)  arc  shown 
in  Figs.  8  and  9.  The  material  constants  for  the  aluminum  were 
identical  to  those  used  above. 

In  Fig.  8  arc  shov;n  deformed  material  coordinates  at 
various  stages  during  the  motion,  vdiilc  in  Fig.  9  arc  shown  cor¬ 
responding  isometric  (Lagrangian)  plots  of  the  hydrostatic  pres¬ 
sure.  The  effects  of  lateral  release  waves  are  clearly  observ¬ 
able. 
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Appendix  A 


One  Dimensional  Finite  Difference  Equations 


A.l  Equations  of  Motion 

Consider  only  discrete  points  in  space  finite 
distances  ^  x.  apart  and  denote  the  initial  coordinate  of 
the  j  **  such  point  x®  .  Similarly,  consider  only  discrete 
times  finite  increments  At  apart,  and  denote  the  n  ***  such 
time  t” .  We  are  then  concerned  with  positions  and  accelerations 

of  these  points  at  these  times, 
denoting  by  x  and  a  f  the 

w  J  j 


J'l 


J*-l 


Acl 


position  and  acceleration  of  the 
j  point  at  the  n  ***  time.  Other  quantities,  such  as  pressure, 
density,  stresses  , etc ., are  averaged  over  the  intervals  Ax:  and 
denoted  by  suppose  that  quantities 

vary  so  slowly  that  linear  interpolation  is  justif ied, e.g. , 

and  similarly  in  time 

.  i  ( J';*'  - 

Such  linear  interpolation  is  in  accord  with  the  approximations 

involved  in  the  finite  difference  analogs  to  partial  derivatives 
correct  to  second  order, i.e., 


A. 2 


A. 3 


V  "  _  V.  ’• 
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11 

H 


t 


ft  *■'/%. 


J _ 

i 


A. 4 


Simplifying  the  notation,  the  momentum  equation  will 

be  \^^:itten 


a 


j_  ^  ( <r  -t  ff) 

P  at 


V\7hcre  we  have  written  a  ~  a*  ^  ^  ~  and 
v;hore  body  forces  have  been  temporarily  omitted.  Applying  the 
above  principles,  the  finite  difference  form  is 


»  -  X 


rk  i  -  xf)  ^  fU  i  x;  -  X.:, ) 


A. 6 


The  velocity  and  coordinate  are  given  by 


J 


i-  a’ 


A. 7 


X 

J 


X 

J 


At 


A. 8 


37 


It  is  simplest  to  take  the  continuity  equation  in  integral 
form,  whence 


where 


A. 9 


is  a  constant  for  each  mesh  point. 

A. 2  Artificial  Viscosity 

In  order  to  link  the  width  of  shock  r.oncs  to  the  mesh 
size,  the  artificial  viscositv  coefficients  are  set  proportional 


to  the 

mesh  size, 

so  tliat  the 

artificial  viscosity  becomes 

=  a:::  I 

JC.  -  X  .  1 

\  J*l  J  / 

"  /i) 

( pJ 

)\f/j  J 

for 

if)>0 

n+i 

=  0 

for  1 

'SI  ^ 

where 

f 

-  a:.j 

1 

f  ■ 

lTic-  sound  speed  docs  not  change  very  rapidly  so  that  the  use  of 

this  equation  docs  not  introduce  any  difficulties. 
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A.  3  Fluid  ConstlLutivc  Equatiion 


Several  equation  of  state  options  are  available.  It 
is  convenient  to  deal  with  a  fljid  separately.  If  the  flow  is 
isentropic,  an  equation  in  the  fomi 

fitled  to  the  isentrope  suffices.  For  more  general  motions,  the 
relation 


cF  =  <F>  *  ^  iF) 

where  /,  and  will  be  expressed  in  terms  of  7 

/,  =  k.  7  /  '  -  *<.  ?  } 

A  =  k,  [  /  +  >.,7  *  k.7*  ^  ■  ••  •] 


A. 12 


is  combined  with  the  energy  equation,  which  in  finite  difference 
form  is 


fj 


n4-i 


A. 13 


where 


Heat  conduction  and  energy  source  terms  have  been  temporarily 
omitted.  Substituting  A. 12  at  t"*’  into  A. 13  provides  an  ex¬ 
plicit  equation  for  the  internal  energy 


I  _  If"**  ^ 

I  i  Jl 


A.  14 
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The  pressure  can  then  be  found  from  equation  A. 12  and  cr *  « ^ 
(P  ^  0  .  The  sound  speed  is  given  by 


I  (/ 


where 


f  0  ^  K  ( ‘'7J{  ^  ^  7  ^  2^^  '  } 


A.  15 


/a/ 


0  -  KI  ^ i 


A, 4  Elastic  Perfectly  Plastic  Constitutive  Equation 


The  stretching  deviator  is 

-  »rV 
-  (*r  - 


+  -t-  ^ 


and  the  elastic  stress  deviator  becomes 


A.  16 


-f'  -  (w)"  * 

a*"  V«  /jf\ 


A. 17 


where  we  will  take  the  shear  modulus  to  depend  on  compression 

*  f-.  f  '  ^  7  *  ?"»  7  ‘  ■  i 


For  ot  *  I  or  3,  the  yield  condition  reduces  to 
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A. 13 


B  = 


Then 

If 

n 

2. 

«  3 

diK 

If 

ff 

>  ^ 

VP 

-  .  §  y 

Then 

0 

^  -  1 

9^  ~  X 

{it ') 

and  the  deviator  stress  work  is  simply 


A. 19 


A. 20 


This  is  used  in  computing  the  pressure  below.  For  o<*  2  a  more 
elaborate  procedure  is  required.  We  compute  the  z  component  of 
stretching  and  stress  deviators 


1  A 

3  / 


A. 21 


4 

e 


2 


n-t-'/x 


G"'"' 

J*<fx 


n  t'/i 


A. 22 


These  equations  are,  of  course  only  valid  if  the  material  is 
elastic,  l^hen  the  material  is  plastic  it  is  necessary  to  solve 
the  set  of  simultaneous  equations  for  the  unknown  stress  compo¬ 
nents  and  proportionality  factor  /  .  Due  to  the  quadratic  form 

of  the  yield  condition,  this  cannot  be  done  explicitly.  Instead, 
a  forward  differencing  scheme  used  by  Wilkins  and  Giroux,"^  and 
Maenchen  and  Sack^is  used.  Referring  to  the  77 plane  representation, 
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v;e  con.sitier  an  elastic  stress 


5  •  2  [  itt‘J 

Then  il 

3  ^  r  y’ 


state  1  at  t",  and  use  equations 
A, 16,  A. 17,  A. 21  and  A. 22  to 
compute  a  new  elastic  stress 
state  2  at  t  .  The  second 
moment  is  then  computed  using 

.  it' .  ii‘  » (U'j'j 


Use 


A  ^ 


A. 23 


since  the  material  is  still  elastic.  However, 

if  f  >  j  y*- 

the  material  has  become  plastic  as  shox/n  in  the  diagram.  Ap¬ 
proximately  the  correct  yield  stress  is  achieved  by  using 


x-zhich  in  elfcct  computes  the  stress  state  at  point  3  on  the 
yield  surface  on  a  radius  vector  from  state  2.  The  procedure 
is  only  valid  x^^hen  the  change  in  Lode's  angle  is  small.  It 
might  be  noted  that  xdicn  «  -  /  or  3 ,  the  symmetry  condition 
t  -  limits  the  attainable  stress  states  to  the 

straight  line  AA'  in  the  diagram  above,  and  the  procedure  given 
for  finding  the  stress  deviators  for  this  case  does  not  involve 

j  f  di  X  z 

this  approximation.  Then  Cp-Z.fl  *  gC  and 
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At 


/I  #'4. 


z  )  "♦'ij 
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A.  25 


Finally  the  hydrostatic  pressure  is  found  as  for  a 
fluid,  but  including  the  deviator  stress  work  in  the  energy 
equation  A. 14 


H  tit  n  -  *»  \ 

”  .  I  ^  0*t  Art 


\  ^  A  S 


.  1  r  "*' 


A.  26 


4±x 

The  pressure  is  found  from  equation  A. 12  and  O'  * 

The  elastic  sound  speed  may  be  taken  as  /rr  times  the  value 
given  in  equation  A. 15 


A. 5  High  Explosive  Constitutive  Equation 

In  order  to  force  the  detonation  front  to  move  at  the 
proper  velocity  the  following  scheme  is  adopted.  The  time  at 
which  burning  is  initiated  in  a  particular  mesh  will  be 


J*</x 


A. 27 


where  x*  is  the  initial  position  of  the  detonation  point  (the 
detonation  is  considered  to  be  initiated  to  the  left,  or  smallest 
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Xj  ),  and  O  is  che  detonation  wave  velocity.  In  orde^  to 
spread  the  detonation  front  over  several  meshes,  a  burn  fraction 
is  defined  as 


=  0 


D 

Ax^ 


( t 


fer  t"'  i  t  ‘ 

A.  28 


I 


or 


where  0  ^  F  ^  I  ,  and  Bj.  is  a  constant  which  determines  the 
thickness  of  the  shock  front.  If  we  now  take 


e 


i- 


A. 29 


the  pressure  is  maintained  zero  until  the  precomputed  detonation 
time,  and  rises  smoothly  until  F=  1.  Since  F  -  I  for  all  times 
thereafter,  equation  A. 29  provides  the  correct  equation  of  state 
of  the  explosion  products  for  the  subsequent  motion. 


Solving  equation  A. 29  together  with  the  energy  equation 
A. 13  for  a  fluid  provides  the  equation  for  the  internal  energy 


L  / F**"'  L ***' 


0*1  0 


I  J.  fz  Ap 


Ap 

~pi. 


A.  30 


The  pressure  is  then  found  from  equation  A, 29. 

Choosing  a  polytropic  law  to  describe  isentropes  and 
using  the  Mie-Grueneisen  assumption  to  generalise  to  other  nearby 
states,  /,  and  can  be  expressed  as 
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A. 31 


K  -  rp 


where  ,  p  and  "y  are  material  constants.  For  the  simple  poly¬ 
tropic  gas,  the  expression  for  the  sound  speed,  equation  A. 15, 
reduces  simply  to 


r 


A.  32 


Other  more  elaborate  forms  for  and  may,  of  course,  be  used 
if  deemed  necessary. 


A. 6  Energy  Checks 

f 

It  is  often  desirable  to  study  momentum  and  energy  of 

*  % 

individual  meshes,  or  summed  over  portions  or  the  whole  of  the 
material  in  the  problem.  The  finite  difference  expressions  be¬ 
low  are  collected  for  convenience,  and  are  incorporated  in  a 
special  subroutine  only  when  required  for  diagnostic  purposes. 


where 


The  mass  in  a  mesh  is  given  by 


M 


k' 


m 


A'- 

1 

for 

oC 

-  1 

k'  = 

ir 

for 

oC 

=  X 

k'  - 

for 

OC 

•t  3 

Thus  the  momentum  in  a  mesh  is 


} 

X  * 


A. 33 


A. 34 
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The  kinetic  energy  is 


K 


•'ji 


J**h. 


A.  35 


while  the  internal  energy  is 


J. 


A.  36 


It  is  useful  to  sum  kinetic  and  internal  energy  over  the  mate¬ 


rial.  Then  if  energy  is  conserved 


A 


Z  J  K  '? 

^*'4  /  j 

£ 


A- '4 
J'f-f/x. 


) 


A,37 


where  is  the  surface  work  done  from  time  t"  to  time 

t"**  (for  example  hy  a  surface  pressure)  and  is  the 

nonmechanical  energy  addition  in  this  time  interval  (for  example, 
chemical  energy  in  an  explosive).  This  provides  a  ver>  useful 
check  on  the  calculation. 


It  is  sometimes  useful,  in  studying  a  particular  motion 
in  detail,  to  consider  the  energy  distribution  in  various  modes. 
Considering  the  time  interval  to  in  each  case,  the 

spherical  elastic  stress  work  done  in  a  mesh  is 


p ^ 


X 


M.  ( 


A.  38 


while  the  spherical  viscous  stress  work  is 


P*p*. 

J*>u. 


A.  39 
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The  plastic  stress  v;ork  can  be  estimated  in  a  ihanner 


consistent  with  the  forxi/ard  differencing  scheme  used  in  .the  con¬ 


stitutive  equation.  Referring  to  the  diagram  in  Section  A. 4, 
note  that  the  plastic  stretching  deviator  is  related  to  the  in¬ 
crement  in  stress  between  states  2  and  3,  i.e., 


‘‘d’ 

P 


2  G 


At 


A.40 


The  plastic  stress  work  per  unit  mass  is  thus  approximately 


p  3^ 


XG  p 


jiy-r)r  A.41 


where  summation  is  implied  over  the  index  a  ,  Using  equation 
A. 24  this  becomes 


I 

2  G  JO  At 


A. 42 


Since  state  3  is  at  yield,  the  second  moment  of  the  stress  tensor 
at  3  is  j  y  •  In  finite  difference  form  the  plastic  work 
done  is  therefore 


P  J*>lx 


S 


s 


11^] 

E  / 


A. 43 


The  total  deviator  stress  work  is 


A^E 


A.  44 
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and  the  elastic  deviator  stress  work  is  simply  the  difference 
between  equations  A. 44  and  A. 43. 

A. 7  Boundary  Conditions 

Boundaries  and  interfaces  require  no  special  treatment 
other  than  to  supply  the  appropriate  quantities  on  either  side  of 
the  boundary  or  interface  in  the  momentum  equation  A. 6. 

For  a  fixed  surface  quantities  <r,  ^ ,  jO  and  in  the 
mesh  outside  the  fixed  surface  are  set  equal  to  the  corresponding 
quantities  iti  the  mesh  inside  the  boundary.  Furthermore  we  set 


where  J  is  the  boundary  index,  solving  for  the  appropriate  x.'* 
outside  the  boundary.  The  acceleration  and  velocity  are  then 
zero  and  is  computed  to  be  equal  to  x*  ,  as  required. 

For  a  free  surface,  quantities  <r,  ^,  />  and  <P  in  the 
mesh  outside  the  free  surface  are  set  equal  to  zero  while 
outside  the  boundary  is  found  as  above.  If  a  surface  pressure 
is  to  be  introduced,  set  <r  =  where  ip  is  sup¬ 

plied  as  input  either  as  an  analytic  fit  or  in  tabular  form. 

For  an  internal  interface  between  different  materials, 
no  special  provisions  are  required  while  the  materials  on  either 
side  of  the  interface  are  in  contact.  However,  when  the  stress 
at  the  interface,  )  where  T  is  the  interface 

index),  exceeds  the  tensile  fracture  stress  of  the  bonding,  two 
free  surfaces  are  formed.  It  is  then  necessary  to  store 
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additional  values  of  x  and  u  for  the  second  surface  and  aibply  the 
free  surface  conditions  when  computing  the  separate  acceljerations 
of  each  surface.  If  at  a  later  time  the  surfaces  attempt*  to  cross, 
(i.e.  X  >  X  where  J-  and  JV  refer  to  the 

left  and  right  hand  free  surfaces  respectively)  the  positions 
and  velocities  of  the  two  surfaces  are  set  equal  and  the  inter¬ 
face  is  treated  as  unseparated  once  more.  Subsequent  separation 
of  the  interface  will  then  occur  at  zero  stress.  Actually  a 
small  (nonzero)  separation  stress  is  used  to  prevent  undesirable 
separation  and  contact  when  the  interface  stress  oscillates 
about  zero.  Fractures  at  the  interior  of  a  material  are  treated 
in  precisely  the  same  way. 

A. 8  Stability  Criterion 

*»♦  V*. 

The  stability  criterion  is  used  to  compute  At 
used  to  advance  ^■he  calculation  on  the  following  time  cycle.  The 
expression  is  evaluated  for  each  mesh  and  the  minimum  is  used  for 
advancing  the  calculation.  Using  a  backward  time  centering,  the 
stability  criterion  becomes 


In  this  expression  it  is  tacitly  assumed  that  the  velocities  at 
j  and  j+i  will  remain  the  same  on  the  next  cycle. 

The  expression,  equation  A. 45,  has  been  found  to  be 
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adequate  for  most  purposes.  Experience  has  shown  that  even  if 
large  velocity  discontinuities  are  introduced  initially ,  as  for 
example  at  the  interface  between  colliding  materials,  the  com¬ 
putation  remains  stable  and  the  discontinuity  is  smoothed  to  the 
normal  shock  width  of  3  or  4  meshes  within  3  or  4  cycles.  How¬ 
ever,  if  an  initial  pressure  discontinuity  is  introduced,  the 
calculation  becomes  unstable  if  equation  A. 45  is  used.  A  heuris¬ 
tic  argument  suggests  the  following  explanation.  A  velocity  dis¬ 
continuity  acts  immediately  to  strengthen  the  stability  criterion 
through  the  velocity  gradient  term  in  the  denominator  of  equation 
Ao43.  A  pressure  discontinuity,  however,  does  not  affect  the 
stability  criterion  until  a  later  cycle  when  the  pressure  dis¬ 
continuity  has  accelerated  the  mesh  points  concerned.  A  scheme 

IX 

to  overcome  this  limitation  has  been  suggested. 


Instead  of  computing  the  stability  criterion  at  the  con¬ 
clusion  of  the  calculations  for  a  particular  cycle,  the  stability 
criterion  is  computed  just  after  computing  the  acceleration 
(equation  A. 6)  but  before  computing  velocity  (equation  A. 7).  A 
forward  time  centering  is  used.  Moving  the  entire  expression 


backwards  one  step  in  time,  we  write 


(it-iBjc'  .  48/  )Au«"'>| 


A. 46 


where  ^  »  etc.,  and  the  spatial  index  is  hence¬ 

forth  omitted.  The  essential  difference  between  equations  A. 45 


and  A. 46  is  in  the  centering  of  the  and  terms.  We  may 
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anticipate  that  and  Au"*’**"  v;ill  differ  naterially 

from  A^"  and  ^u"~'^if  the  acceleration  gradient  is  high. 
Such  an  acceleration  gradient  would  result  from  a  large  pres¬ 
sure  discontinuity.  The  sonic  velocity  c  by  comparison  does 
not  change  drastically  even  with  large  changes  in  pressure  so 
that  it  is  sufficient  to  use  c " .  At  this  point  in  the  calcu¬ 
lation  the  velocities  and  position  are  not  yet  available  at 
n+l  and  '/x  respectively.  However,  a  good  estimate  is 


u  =  U 

J  J 


X  V 

J 


«  AC 


A. 47 


j  '  j 

since  At  will  not  be  expected  to  vary  too  drastically  from 


cycle  to  cycle, 


Inserting  equations  A. 47  into  A. 46  and  investigating 
limiti^.g  solutions  for  large  and  small  values  of  Ao^”  leads  to 
the  sufficient  condition 

^  A./^8 

i-l  )Ia  U*'*'*!  +  +  Ak"  IAh’'!  j 

where  **»  x  -x."  ,  etc. 

A  constant  K  is  inserted  as  a  convenience  for  strength¬ 


ening  the  stability  criterion  artificially  if  this  should  be 
desirable . 


Appendix  B 


Two-Dimensional  Finite  Difference  Equations 

B.l  Equations  of  Motion 

Consider  discrete  points  in  space  formed  by  the 

inLcrsccLicn  oT  maLerial  coordinate  lines  distances  A>i  and 
A  2  apart.  Similarly,  consider  only  discrete  times  finite 
increments  AX  apart.  Wc  use  three  indices  to  denote  values  of 

of  the  2-  material  coordinate 
at  the  n*  time,  e.g.  <3.  i  : 
Quantities  such  as  pressure, 
density,  stress,  etc.  are  av¬ 
eraged  over  the  meshes  formed 
by  the  finite  difference  grid, 
and  arc  denoted  by  p.  " 
etc . 

indexed  in  this  way  in  the  com¬ 
puter,  in  writing  down  the  finite  difference  equations  it  is 
more  convenient  to  use  the  notation  shown  in  the  sketch. 

Suppose  that  quantities  vary  so  slowly  in  time  that 
linear  interpolation  is  justified  ,i .e. , 

=  a  c  ’ 

and  similarly  in  space.  Such  linear  interpolation  is  in  accord 


quantities  at  the  intersection 
and  X  maLerial  coordinate 


VJliile  quantities  are 


52 


v/ith  the  approximations  involved  in  the  finite  difference  ana¬ 
logs  to  partial  derivatives,  correct  to  second  order,  i.e.  for 
time  derivatives  at  0 


ll 
a  i 
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n  ♦ 
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n  4 


B.2 


Finite  difference  analogs  to  spatial  derivative  are  obtained  by 
the  use  of  Green's  Transformation.  For  quantities  which  are  av¬ 
eraged  over  the  meshes,  Green's  Transformation  applied  to  the 
circuit  A  B  C  D  gives  the  gradients  at  0 


I  ii* 
fa* 


;  *  IK 
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±  if 
f  ^2 


where  A  is  the  area  of  the  quadrilateral  A  B  C  D. 

The  momentum  equations,  temporarily  omitting  body 
forces,  become 
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where  we  use  equations  B,3  to  represent  the  derivative  terms 
and  the  equations  are  centered  at  0  in  space  at  t"  in  time. 

The  term  appearing  in  these  equations  is  required  at  0, 

but  p  and  fl  are  quantities  averaged  over  the  meshes.  The  sim¬ 
plest  solution  is  to  use 

If^).  ■  S  A  +  A  *  A  ^*)  B-5 

For  cases  where  areas  of  severe  distortion  must  be  handled  as 
realistically  as  possible,  it  is  better  to  use  the  more  accurate 
but  much  lengthier  relation 

-  I  P,  *  A  *  Ps  *  A 

B.6 


is  the  area  of  triangle  A  0  D  etc.  We  have  used  equation  B.5. 


The  last  terms  on  the  right  of  equations  B.4  also  re¬ 
quire  interpolation.  A  variety  of  interpolation  schemes  are 
possible,  we  have  used 


1 

p>c 

where  /? 
later. 


=  i  ^  ^  II'.  ^  J  B.7 

and  m  are  the  areas  and  masses  of  the  meshes,  defined 
This  allows  a  relatively  simple  treatment  at  boundaries. 
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The  velocities  and  coordinates  of  cell  vertices,  using 


equation  B.l,  are  given  by 
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B.2  Continuity  Equation 

The  continuit}  equation  is  taken  in  integral  form.  The 
volume  of  mesh  3  is  given  within  a  factor  '  by 


where  and  are  the  areas  of  the  triangles  B  F  C  and 

B  0  C  respectively  at  time  i  " 

At_  =  I  f  (^c  -  ~  ^f)  -  -  ZftJ J 

B.ll 

A  u  ~  Z  j  [  ^0  “  ^o)  ~  ^5  "  ^o)  J 
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and  and  are  the  centroids  of  the  triangles  B  F  C  and 

B  0  C  respectively  at  time  t  ” 
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The  continuity  equation  is 


B.13 

where  is  a  constant  for  each  mesh,  evaluated  at  time 

m,  =  jo/  I  fiL  B.14 


nicv  area  of  mesh  3  is  simply 


B.15 


B.3  Stretching  and  Spin 


The  components  of  the  stretching  deviator  and  spin 
tensor  are  given  by 
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InLcrpolation  in  botb  space  and  time  is  required  in  the  finite 


difference  analogs  to  the  velocity  gradient  terms.  Applying 
Green's  Theorem  to  the  circuit  B  F  C  0  and  using  linear  inter¬ 
polation,  we  get 


y'  ^  _ 


^  fit-. 


ij _ -I 


^2 


-z; 


which  is  centered  at  3  and  time  t  .  Also 
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B.4  Artificial  Viscosity 

In  order  to  link  the  width  of  shock  zones  to  the  mesh 
size,  the  artificial  viscosity  coefficients  arc  set  proportional 
to  the  mesh  size.  The  bulk  viscosity  becomes,  with  all  quantities 
centered  at  3 
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A 


I  -  [b,  /r^  (fj]  ‘  } 

/  B.19 

for  y>  ^  ^ 

<^**'-  0  for  ^  <  0 

The  cicviator  viscous  stresses  take  the  form 

i\’T'  ’  f"''l '*•20 

vvhere 

*  '  Ma«  (  X,  ,  X,,  fx,  ,  X,  ,x,, 

z  '  May(z.,  z,,  z^,  -  t1in.(z,, 
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and  arc  scL  equal  to  zero  when  the  corresponding 

stretching  dcv^iators  arc  positive.  These  expressions  are  not 
properly  centered,  Init  s.ince  the  density  anc'  sonic  velocity  do 
not  change  rapidl}-,  tliis  does  not  introduce  any  dif  I  icul  tics . 

The  v/orlc  done  hy  the  viscous  stresses  in  the  time  in¬ 
terval  At''*  is  given  by 


where 


f 


(f 


n*i 


for  the  spherical  part,  and 


B.21 


for  the  deviator  part. 
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B.5  Stress  Deviator 


VJhcn  the  material  is  elastic,  the  stress  deviators 
are  found  from  the  stretching  deviators  by  the  differential  equa¬ 
tions 


1  u; 
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at  ^  ^ 
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B.23 


Tlie  differential  terms  on  the  left  arc  set  into  finite  differccc 
form  using  equation  B.2,  while  the  other  stress  terms  on  the  left 
must  be  interpolated  to  t using  equation  B.l.  This  leads 
to  the  implicit  set  of  equations 


iUT  -irf 
irr-(xr^ 


-  At  -  2  At  G  -d” 
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where  all  quantities  except  stresses  are  centered  at  t 

If  we  choose  the  time  increment  At 
small  so  that 


(Ai"* 


r  «/ 


sufficiently 


B.25 


then  these  equations  may  be  expressed  in  a  relatively  simple  ex¬ 
plicit  form.  The  condition  B.25  implies  that  the  rotation  from 
t "  to  ^  ig  small.  Since  the  spin  is  the  angular  velo¬ 
city,  the  rotation  in  time  A  t  ig 

w'*  B.26 

Defining 


j.  4  4,  XX 


3.  At"”"'  C,  {‘‘d") 


n  ♦ 


A:i 


2  At  G  '*"*  'V 


B.27 


*  3  At"”"” 

where  we  take  the  shear  modulus  as  a  function  of  compression 

^  *  i-l  '  *  f'  1  *  r-  7'  *  ■■■  }  B-28 

we  find  that  equations  B.24  solved  for  the  stress  deviators  at 
time  t  become  approximately 
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It  is  readily  seen  that  in  each  equation  the  second 
term  on  the  right  represents  the  stress  increment  due  to  stretch 
ing,  while  the  last  term  represents  that  due  to  a  small  rotation 

'fliesc  equations  are,  of  course,  only  valid  when  the 
material  remains  elastic.  When  the  material  is  plastic,  only 
the  clastic  stretching  should  appear  in  the  above  equations, 
and  it  is  necessary  to  solve  the  set  of  simultaneous  equations 
for  the  unknov/n  stress  components  and  proportionality  factor 
A  method  of  forward  differencing  used  by  Wilkins  and  Giroux,^ and 
Maenchen  and  Sack^is  used.  Consider  an  elastic  stress  state  at 
time  .  If  we  v/ere  to  transform  to  principal  coordinates, 
this  stress  state  could  be  plotted  in  the  IT  plane  representa¬ 
tion  as  point  1,  say  (see  diagram  on  next  page). 
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Equations  B.2o  to  B.2o  arc  used 
to  compute  a  nev;  elastic  stress 
state  at  ^  ,  Transforming 

to  a  nev;  set  of  principal  co¬ 
ordinates  (which  arc  in  general 
rotated  with  respect  to  the 
previous  ones)  v;e  can  again 
represent  this  state  in  the  IT 


plane.  The  second  moment  is 


B.30 


Then  if 


3  *  jy" 


Use 


B.3i 


since  the  material  is  still  elastic. 

Hov;ever  if  ]J  >  y  ^  ^ 

the  material  has  become  plastic  as  shown  in  the  diagram.  Ap¬ 
proximately  the  correct  yield  stress  is  achieved  by  using 
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B.32 


which  in  effect  computes  the  stress  state  at  point  3  on  the  yield 
surface  on  a  radius  vector  from  state  2.  The  proceedure  is  only 
valid  v;hon  the  change  in  Lode’s  angle  is  small.  Due  to  the  ro¬ 
tation  of  the  principal  axes,  the  principal  stress  components 
must  be  found  before  Lode's  angle  may  be  determined.  We  have 
for  principal  stress  deviators 


{y  * i /(u"  - ^ 


B.33 


where  a  *  /,  2  respectively  v/hen  the  upper  or  lower  sign  is 
used.  Tlie  angles  of  the  principal  axes  with  respect  to  the  K 
coordinate  axis  are 
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a  rc 


tat 
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and  Lode's  angle  is  given  by 


taf* 


/T 


it'  -  it" 

U‘  * 


B.35 


Note  that  principal  components  and  directions  are  not 
required  except  optionally  for  diagnostic  purposes. 

Finally  the  deviator  stress  work  is  given  by 


^  ^  f  I 


[(XT'*(irfj[ 


B.6  Pressure 

The  pressure  is  taken  as  a  function  of  the  thermody¬ 
namic  state  in  the  form 

eP  ’  ilP)  *  ^  AOV  B.37 

where  in  terms  of  7  ‘  (f 

A,7  .  7* .  j 

L  •  {  I  *  h,>i  *  A.7^  *  •••  j 
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v/hilc  the  energy  equation  in  iinite  diflerence  form  is 


r'  .r  (,r' 

where  4%o*  Is  given  follov;ing  equation  B.21.  Solving  equation 
B.37  and  Bo38  for  <5'"’' 


The  pressure  can  then  be  found  from  equation  B.37.  Then  the 
total  stress  components  are  given  by 

t"  =  ti’"  i-  V”  -  *  t) 

-  tr  *  V“  *r) 

B.40 

t  -  J  *  f 


f“  -  t"  =  I'f"  i-  tt”  *  Z  ■‘f’"  *  *f^"- 

Tlie  bulk  modulus  K  is  obtained  by  differentiating  equation  B.37 
v;ith  respect  to  J>  at  constant  entropy.  Assuming  that  Poisson's 
ratio  remains  near  1/3,  a  good  expression  for  the  elastic  sound 
speed  is  C  .  /is  K/j>  ■  whence 
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where 


f  f  '  *  i**.  7  *  -5  }  (I  -fj 

f  -  /)„[  A,  ■»■  iA,7  +  JAj  ?*+  •  •■ }(/-?; 


B.7  Energy  Checks 


It  is  often  desirable  to  study  momentum  and  energy  of 
individual  meshes,  or  summed  over  portions  or  the  whole  of  the 
material  in  the  problem.  The  finite  difference  expressions  be¬ 
low  are  collected  for  convenience,  and  are  incorporated  in  a 
special  subroutine  only  when  required  for  diagnostic  purposes. 


The  mass  of  mesh  3  is  given  by 

0(~l 


M,  =  (Xrr)  '  mj 


Thus  the  momentum  components  in  mesh  3  are  approximately 
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M. 
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U 


3  a  -  3 

while  the  kinetic  energy  components  are 


K«  "/'*  .  i  (  u«  " 


B.42 


B.43 


B.44 
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where  in  equations  B.43  and  B.44  we  interpolate  the  velocity  com¬ 
ponents  as 


f  *  U,  *  J 


B.45 


The  internal,  energy  is 


£ 


*♦*4  . 
J 


B.46 


It  is  useful  to  sum  kinetic  and  internal  energy  over  the  mate 
rial.  Then  if  energy  is  conserved 


sJ 


k* 

i,J 


n-t  '4 


B.47 


where  is  the  surface  work  done  from  time  t "  to 

time  t  (for  example  by  a  surface  pressure)  and 

is  the  nonmechanical  energy  addition  in  this  time  interval  (for 
example,  chemical  energy  in  an  explosive).  This  provides  a  very 
useful  check  on  the  calculation. 


It  is  very  simple  to  find  the  energy  dissipation  in 
various  modes  during  the  time  interval  to  f  ,  Xhe  spher 

leal  elastic  stress  work  is 


3 


B.48 


while  the  spherical  and  deviator  viscous  stress  work  are  given  by 
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B.49 


3 


^3 


a 


The  plastic  stress  work  is,  as  for  the  one  dimensional  case, 

T  J 


3 


U  Mj 


ri 

J  ^  _ 

IT 


ar^(pr  *■/>:) 

The  total  deviator  stress  work  is 


B.50 


A-e 


B.51 


and  the  elastic  deviator  stress  work  is  simply  the  difference 
between  equations  B,  51  and  B.50. 


B.8  Boundary  Conditions 

Boundaries  and  interfaces  require  no  special  treatment 
other  than  to  supply  the  appropriate  quantities  in  the  meshes 

surrounding  a  boundary  or  interface  point  in  the  mo^hentum  equa- 

/ 

tions  B.4. 

For  a  fixed  surface  along  which  the  material  may  slide, 
normal  stresses  t  ***  and  t  ,  densities  areas  ind  masses 
h  and  m  in  the  meshes  outside  the  boundary  are  set  equal  to  the 
values  in  the  corresponding  meshes  inside  the  boundary,  while 
the  shear  stresses  t  **  are  reflected  antisymmetrically ^  It  is 
also  necessary  to  supply  the  coordinates  of  mesh  vertices  outside 
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the  boundary.  If  A  C  is  the 
boundary,  we  require  coordinates 
of  point  D  outside  the  boundary 
which  is  a  reflection  of  point 
B  inside  the  boundary.  The 
equations  are 

*2^} 

B.52 

where 


—  /  Z—  /. 

;  A*  9  I 


^ra  ‘ 


•  ^  *  ^/n  * 

For  a  free  surface,  the  stresses  and  densities  are  set 

equal  to  zero  in  the  meshes  outside  the  boundary.  To  prevent 
difficulties  with  zeros  in  the  denominators  of  terms  in  equation 
B.7  , ^ and Z for  the  cells  outside  the  boundary  are  treated 
as  for  a  fixed  surface.  If  a  surface  pressure  is  to  be  intro¬ 
duced,  set  t“"*  t**»  ~~  iP  where  (tj  is  sup¬ 

plied  as  input  either  as  an  analytic  fit  or  in  tabular  form. 


Corners  require  no  special  consideration,  except  that 


n 

1 - 1 


I 


side  the  boundaries  (mesh  4  in 


,  7L„  are  used  in  place  of 
Xj,  2g  to  find  the  position 
of  point  D  outside  the  boundary 
(and  similarly  for  point  C) . 
Values  in  the  comer  mesh  out- 
the  diagram)  may  be  obtained 
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either  by  appl-ying  the  proper  boundary  condition  along  A  0  to 
find  values  in  rneSh  1,  and  subsequently  applying  the  proper 
boundary  condition  along  B  D  to  find  values  in  meshes  3  and  4, 
or  by  applying  the  proper  boundary  conditions  along  B  0  to  find 
values  in  mesh  3,  subsequently  applying  the  proper  boundary  con¬ 
ditions  along  A  C  to  find  values  in  meshes  1  and  4.  In  either 
case  identical  results  are  obtained.  Reentrant  corners  are  gen- 


B 


erally  formed  by  the  meeting  of 
two  free  surfaces,  and  it  is  un¬ 
necessary  to  find  values  of  x. 
and  2  at  points  C  and  D  (in  the 
c  diagram)  by  equation  B<.S2, 

However  no  harm  is  done  if  equation  B.52  is  used,  and  it  is  more 
convenient  to  retain  complete  generality  at  all  boundary  points. 


Q 

•3 

1 

1 

•  “f  1 
1 

- 1 

Interfaces  between  different  materials  require  no  spe¬ 
cial  provisions  while  the  materials  remain  in  contact  without 
sliding.  Sliding  interfaces  and  separated  interfaces  require 
individual  treatment,  as  do  internal  fractures.  These  can  be 
incorporated  as  required. 


B.9  Stability  Criterion 

The  stability  criterion  is  used  to  compute  ^ 

which  is  used  to  advance  the  calculation  on  the  following  time 
cycle.  Using  a  backward  time  centering 
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K //)*"' 


B.53 


where  the  minimum  is  taken  over  all  meshes o  A  constant  K  is 
inserted  as  a  convenience  for  strengthening  the  stability  cri¬ 
terion  artificially  if  this  should  be  desirable. 
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Appendix  C 


Evaluation  of  Constitutive  Equation  Constants 

The  constitutive  equations  appearing  in  Section  III 
have  been  specialised  to  specific  classes  of  materials,  by  mak¬ 
ing  assumptions  that  the  materials  are  perfect  fluids  or  com¬ 
pressible  elastic-perfectly  plastic  solids.  The  assumption  was 
made  (equation  3.6)  that  the  spherical  part  took  the  form 

r>  -  fjfii  *  £.  L(j>)  c.i 

where  and  were  left  as  arbitrary  functions,  although  it 
was  intimated  that  power  series  expansions  could  be  used.  It 
is  not  necessary  to  limit  ourselves  to  the  form  of  equation  C.I, 
and  it  is  quite  feasible  to  use  a  more  general  expression 

P  -/(/>.<?;  C.2 

although  the  convenient  explicit  finite  difference  scheme  for 
finding  the  internal  Cinergy  and  pressure  outlined  in  Apendices 
A  and  B,  which  is  made  possible  by  the  form  of  equation  C.l^must 
be  abandoned,  and  replaced  by  an  iteration  scheme.  It  was  found 
that  an  iteration  scheme  often  required  more  than  three  itera¬ 
tions  to  assure  accuracy,  ^nd  it  has  been  deemed  preferable  to 
express  the  spherical  relation  in  the  form  of  equation  C.I  where- 
ever  possible. 

The  present  appendix  is  concerned  with  evaluating  some 
special  forms  of  the  functions  and  /^for  some  restricted  cases, 
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and  is  thus  much  more  specialised  and  restricted  than  the 
material  in  the  body  of  the  report.  The  functions  are  developed 
for  a  solid  or  fluid  from  measured  Hugoniot  data  using  the 
Mie-Grueneisen  equation.  For  a  gas,  the  functions  are  developed 
from  measured  Chapman- Jouguet  isentropes  by  somewhat  similar 
means.  These  approaches  should  be  reasonably  accurate  up  to 
moderately  high  pressures,  such  as  those  encountered  in 
solid-explosive  systems.  Other  means  of  evaluating  the  spheri¬ 
cal  constitutive  equation  could,  of  course,  be  substituted. 

C.l  Solid  or  Fluid 

Hugoniots  (loci  of  states  attainable  in  a  single  shock 
compression  from  the  normal  state)  have  been  evaluated  for  a 
wide  variety  of  materials  from  experimental  determinations  of 
shock  and  particle  veloci  ■ies*^by  means  of  the  Rankine-Hugoniot 
relations 

[a] 

7  *  W 

W  '  P.  C.3 

LpI? 

i  fic 

where  p  -  )/P  ,  (  refers  to  the  state  ahead  of 

the  shock,  []  refers  to  the  jump  in  variables  across  the 
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/ 


shock,  and  U  is  Lhe  shock  velocity  (  [WJ  •  U  -  u„  ). 
Results  are  often  fitted  to  the  power  series. 

=  k,  7  (  I  ♦  k.'^  *  l<'  '  *  ■  )  C.4 


where  the  normal  state  is  taken  at  ^  O  • 

For  moderate  compressions  it  is  found  experimentally*'^ 
that  the  shock  velocity  is  a  linear  function  of  particle  velocity 

U  =  c  -h  s  u  C.5 


where  c  is  the  bulk  sound  speed  and  s  is  a  parameter 
terial  ahead  of  the  shock  is  taken  to  be  stationary, 
so  that  equations  C,3  a  and  b  become 


(the  ma- . 

Wo  -  ^  ) 


C.6 


Expanding  to  the  form  of  equation  C.4j  we  obtain  immediately 


2s 


C.7 


Terms  in  kj  and  higher  need  not  be  retained  since  the  linear 
relation,  equation  C.5, is  not  valid  except  when  y  is  small, 
(This  follows  directly  from  equation  C.6,  for  we  note  that 
p  when  -►  l/5  ,  a  physically  unreasonable 

situation.  Since  s  is  generally  found  to  lie  between  1  and  2 
for  most  materials,  equations  C.6  and  C.7  must  be  limited  to 
compressions  of  perhaps  0.1.) 

Note  that  the  coefficient  k'„  represents  the 


7  = 


adiabatic  bulk  modulus  of  the  material  at  zero  pressure. 
The  Murnaghan  equation 


has  occasionally  been  used  to  fit  experimental  shock  Hugoniots/^ 
Expanding  to  the  form  of  equation  C.7  we  obtain 

*  /)/  k.'-  s(f*l)  kj C.9 

where  again  higher  order  terms  need  not  be  retained. 

It  might  be  noted  that  the  assumption  of  a  linear 

Hookean  material  with  Lam^  constants  X  and  n  leads  via  finite 
\  ^ 

strain  theory  to  a  relation  between  the  pressure  and  volumetric 

^  .  \€ 
strain 

C.IO 

where  p  »  /«  (  f  ~ 

which  is  easily  expanded  into  the  form  of  equation  C.7,  whence 

K'  ‘  *  j k/  *  ki  *  C.ll 

The  above  relations  provide  a  direct  means  of  compari¬ 
son  of  experimental  data  fitted  to  the  different  empirical  ex¬ 
pressions,  and  allow  the  data  to  be  easily  put  into  the  form  of 
equation  C.7,  which  is  used  as  a  basis  for  the  following  dis¬ 
cussion  . 
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The  Hugoniot  data  is  generalised  by  making  the  as¬ 
sumption  that  the  pressure  offset  of  a  state  from,  the  Hugoniot 
is  proportional  to  the  energy  density  offset 


[r  '  r***]  *  ^h)  C.12 

where  y  is  the  Gruneisen  ratio.  Ph  (p)  and  (f)  are 

available  from  equations  C,7  and  C.3c,  while  '^  ( )  has  been 
evaluated  approximately  for  numerous  metals  in  the  form'* 


C.13 


The  thermodynamic  relation  at  zero  pressure  should  be  noted 

T(.  •  C.14 

where  <x  is  the  volumetric  thermal  expansion  coefficient  and 
is  the  specific  heat  at  constant  pressure.  For  moderate  com¬ 
pressions,  it  seems  to  be  sufficiently  accurate  to  take  7  * 
constant. 


Equation  C.12  can  easily  be  put  into  the  form  of  equa¬ 
tion  C.l, whence 

/i  '  Ph  -  ^ 

C.15 

K  ^  y  P 

Inserting  equations  C.3  and  C.13  for  and  7  allows 

these  functions  to  be  put  in  the  form  of  power  series  in  7 

k,  7  *  K>1^  *  *  ■■  )  C.16 

K  ( I  ♦  *>.7  *  ♦  K’l**  ■■■  ) 
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where 

.  k/ 

k,  -  k, '  -  i 

K.  -  kJ-iK(K*k.') 

h,  -  k,  -  »  y,  ♦  A,  k,'  ♦  k,'^ 

etc . 


k.  •  />c  y. 

h,  »  1*7, 

K  I  *  X 

k,  •  I  *  7,  *  7^  +7j 


For  a  solid,  the  deviator  part  of  the  constitutive 
equation  depends  on  the  shear  modulus  which  was  also  expressed 
as  a  power  series  in  the  compression.  The  variation  in  shear 
modulus  is  considered  in  another  report.  Since  velocity  of  lon¬ 
gitudinal  elastic  waves  is  given  by 

•  /(X*  i  e>)/f'  C.17 


the  error  in  elastic  wave  velocity  due  to  an  error  in  d  is  only 

or  for  Poisson's  ratio  near  j  ,  ^  ^  ^  •  Also, 

the  error  in  stress  deviator  is  proportional  to  the  error  in  d  , 
but  in  the  presence  of  a  pressure  component,  the  error  in  the 
total  stress  is  smaller  than  this.  For  moderate  compressions, 
it  is  therefore  probably  sufficient  tn  take  d  as  a  constant. 

Jones  has  investigated  the  variation  in  d  on  the  basis  of  several 
assumptions,  and  has  discussed  means  of  finding  higher  order 


terms  in  a  series  expansion  for  d  . 


17 
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C.2  Nonideal  Gas 


The  spherical  constitutive  equation  for  a  gas  must 
take  a  different  form  from  that  for  a  solid  since  in  the  latter 
the  density  is  nonzero  at  zero  pressure,  and  pressures  may  be¬ 
come  negative  (or  tensile). 


We  consider  only  a  gas  in  which  a  known  isentrope  may 
be  fitted  by  a  polytropic  law 


= 


C.19 


where  /.  »  is  a  constant  evaluated  at  a  knowri  point. 

This  approach  is  applicable  to  a  simple  description  of  gaseous 
explosion  productsl^  Chapman-Jouguet  isentropes  have  been  meas¬ 
ured  for  some  explosives,  and  />^  are  then  the  C7  pressure 

and  density  respectively. 


The  isentrope  given  by  equation  C.19  can  be  expanded 
to  cover  nearby  states  by  again  assuming  that  the  pressure  off¬ 
set  is  proportional  to  the  energy  density  offset.’^ 

(p-Pi)  =  7j,  (£  -  Si)  C.20 

In  this  case  is  the  internal  energy  on  the  isentrope,  given 
by 

C.21 

Equation  C.20  can  then  be  put  into  the  form 

C.22 
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where 


/i  ■  ^  f'’ 

k  ‘  yp 


where  we  have  substituted  from  equations  C.19  and  C.21,  and 
where 
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Stress  Profiles  at 


1  Continued 
6  and  8  microseconds 
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Figure  2 


Fluid  Plate  Impact  at  Initial  Velocity  1. 
Stress  Profiles  at  2  and  4  microseconds 


9  km/sec 


Figure  2  Continued 

Stress  Profiles  at  6  and  8  microseconds 
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Figure  2  Concluded 

Stress  Profiles  at  10  and  12  microseconds 
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Figure  3.  Free  Surface  Velocity  as  a  Function  of  Target  Plate  Thickness  in  Terms 
of  Driver  Plate  Thickness  for  a  Driver  Plate  Velocity  of  1.9  km/sec 
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Figure  4.  Free  Surface  Velocity  as  a  Function  of  Target  Plate  Thickness  in  Terms 
of  Driver  Plate  Thickness  for  a  Driver  Plate  Velocity  of  1.2  km/sec 
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Figure  5.  Equation  of  State  Schematic 
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Figure  6./  Stress  Profiles  in  an  Aluminum  Sphere  due  to 
/  a  Spherical  Cavity  Change  of  Pentolite 
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Figure  7.  Cavity  Radius  vs.  Time 
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Figure  8.  Deformed  Mesh  Shapes  for  the  Impact  of  an  Aluminum 
Cylinder  on  a  Smooth  Wall 


Figure  8.  Deformed  Mesh  Shapes  for  the  Impact  of  an  Aluminum 
Cylinder  on  a  Smooth  Wall 
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